#!/usr/bin/env python2.6
# -*- coding: utf-8 -*-
# Copyright (c) 2008 Qtrac Ltd. All rights reserved.
# This program or module is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version. It is provided for educational
# purposes and is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.


import sys
import pickle
import string
import csv
import os
import io
import numpy as np
from scipy import stats
import time
import entropy_shannon
import entropy_shannon_weight

"""Trimming the dataraw file to be used for the entire procedure"""

#15. JANUAR ELLER DEROMKRING 2012
#NEDENFOR HAR JEG PRØVET AT INDICERE HVOR HVILKE SCRIPT SKAL PALCERES


#DET FORLANGES PT AT BRUGER SUPPLERER DE TRE FILES SOM CSV (SENERE MÅSKE ANDRE TYPER EG EXCELL, DAT MM.


#BEGGE DO VIL SAMLET RESULTERE I EN ID-LISTE MED SUBJECTS DER OPFYLDER BEGGE OVENSTÅENDE, DVS FILERNE KAN BLIVE STRÆKT REDUCEREDE. 
#SPØRGSMLET ER: BRUGER LOADER SINE FILER, MEN EVT. (HELT SIKKERT) BRUGE DATA PÅ FLERE FORSKELLIGE MÅDER, FEKS I ET EXPERIMENT KUN BRUGE KØN OG CHOLESTEROL IT ET FORSØG, BRUGE CLASSERE 1, 2 OG 3 I ET ANDET ETC. ALTSÅ, BRUGER SKAL KUNNE LÆGGE SINE RÅFILER IND I EN MAPPE I PROJEKTET OG HEREFTER MANIPULERE SOM ØNSKET; DET KRÆVER VEL EN IMPORT DER HEDDR "RAW-FILES" (INCL GENER) OG HERFRA "CONSTRUCT PROJECT", HVOR OVENNÆVNTE IMPLEMENTERES.


#GENER ELLER SNPS: DER ER TO MULIGHEDER
#	1) RÅ-FILEN KAN STRIPPES FOR ALLE SNPS DER ER MONOMORFE UANSET DO. DET ER SÆRDELES RELEVANT NÅR DER MANGE SNPS. DET SKAL SKE I SERVEREN, MEN HVOR SKAL FILEN PLACERES? JEG TROR I RÅ-MAPPEN MED ET NAVN SOM ASÅ VÆLGES I "CONSTRUCT FILES".
#	2) MED DE ID DER ER TILBAGE EFTER PARTITION-TRAIT PROCEDUREN KAN NOGLE SNPS VÆRE MONOMORFE. DISSE FILTRERES FRA PÅ SERVEREN. 

#VI SKAL SÅLEDES STÅ TILBAGE MED FILTREREDE FILER I "CONSTRUCT PROJECT".


def filter_genotypes(genotype_matrix, genotype_names):  
        """
        filter monomorphic SNPS from entrie snp-file
        """
        #1. call from 
        monomorphic_genotypes = []

        mono_discard = []#list of itmes (snps) to be deleted from matrix

        num_subs = genotype_matrix.shape[1]

        print("no subj") 
        print num_subs

        for snp_index in xrange(len(genotype_matrix)):
            for sub_gt in xrange(num_subs):
                if genotype_matrix[snp_index,0] != genotype_matrix[snp_index, sub_gt] \
                    and genotype_matrix[snp_index, sub_gt] != -1 and genotype_matrix[snp_index, 0] != -1:
                        # not monomorph
                        #print genotype_matrix[snp_index,0], genotype_matrix[snp_index, 0:sub_gt+1]
                        #OK MEN ER JO LIDT MERE LIGETIL I GUIEN
                        break       
                if sub_gt == num_subs-1:
#                    print("mono")
 #                   print sub_gt,num_subs
  					#Just checking that all subjects have been checke ie the end of file is reached
                    mono_discard.append(snp_index)

                    #print "found monomporph", genotype_names[snp_index]
                    if genotype_names:
                        monomorphic_genotypes.append(genotype_names[snp_index])
                        #genotype_names.remove(genotype_names[snp_index])
        
        filtered_matrix = np.delete(genotype_matrix, mono_discard,0)
        mono_discard.reverse()#WHY??to remove frem right?
        for i in mono_discard:
            genotype_names.pop(i)
        return filtered_matrix, genotype_names, monomorphic_genotypes

def filter_genotypes_select(selection_matrix,snps_not_HWE,job_value_temp
):  
    """
    filter monomorphic SNPS from selected cleaned 3-way trait-snp file
    """

    genotypeA = [0,0,0]
    genotypeB = [0,0,0]

    num_subs = selection_matrix.shape[0]

    for gt in xrange(len(selection_matrix)):
        genotypeA[int(selection_matrix[gt,1])-1] +=1
        genotypeB[int(selection_matrix[gt,2])-1] +=1

    gtzeroA = 0
    gtzeroB = 0 
    for snpA in range(0,len(genotypeA)):                     
        if genotypeA[snpA] == 0:
           gtzeroA +=1
    if gtzeroA == 2:
       gene_valid = 1
       return gene_valid
    for snpB in range(0,len(genotypeB)):                     
        if genotypeB[snpB] == 0:
           gtzeroB +=1
    if gtzeroB == 2:
       gene_valid = 2
       snps_not_HWE.append(job_value_temp["snp2"]) 
       return gene_valid


def epi_selection(genotype_names, genotypes, phenotype_names, phenotypes, significance_level, bonf_correct,analysis_output):
    
    results = []
            
    snps_in_file= len(genotype_names)
    
    if snps_in_file < 2:
        return results 
    
    
    num_snp_tests = snps_in_file*(snps_in_file-1)/2 # b: number of snp combinations

    #*SNP og gene name i dict
    job_value_temp = {}
#her eller længere nede?
    job_value_temp["snp1"] = "---"	
    job_value_temp["gene1"] = "---"	
    job_value_temp["snp2"] = "---"	
    job_value_temp["gene2"] = "---"	
    job_value_temp["PhenVar"] = 0
    job_value_temp["TgenVar"] = 0
    job_value_temp["BonFer"] = analysis_output["bonferroni_correction"]
    job_value_temp["BonPerform"] = analysis_output["perform_bonf"] 
    job_value_temp["Significance"]=analysis_output["significance_level"]
    job_value_temp["MAF"]=analysis_output["minimumAllFre"]
    
    if job_value_temp["BonPerform"] == True:
        cut_off = float(job_value_temp["BonFer"])
        #print cut_off    
    else:
        cut_off = float(job_value_temp["Significance"])
       
    for trait_index in range(0,len(phenotype_names)):
        name_trait = phenotype_names[trait_index] 
        print(name_trait)
        #trait_use = []
        epi_sum_out = []
        trait_result = {}
#        entropy_result ={}
        trait_result["trait_name"] = name_trait
        trait_result["snp_interactions"] = []
        trait_result["AllGenFreq"] = []
        trait_result["Entropy"] = [] 
        trait_result["Entropy_data"] = []     
        trait_result["Entropy_data_weight"] = []                        

#NEDENSTÅENDE CALCULATIONER ER CENTRALE OG BEREGNES PÅ SERVEREN, PASSES TIL GRIDDEN
      
        trait_mean = np.average(phenotypes[trait_index])
        trait_var = np.var(phenotypes[trait_index])

        #print(trait_var_class)   
        trait_result["trait_mean"] = trait_mean
        trait_result["trait_variance"] = trait_var
        
        snps_not_HWE = []#in use?


#check first snp here; worthwile ???
        first_snp = genotype_names[0]
#        print first_snp

        for snp1_index in range(0,len(genotype_names)):
#            name_snp1, gene_snp1 = genotype_names[snp1_index]
#####NBNBNB This is just a rescue-operation as call for gene-names somtimes turn up 2 names. Has to be solved in eftech. Do line can be used problem is fixed.
            name_snp1 = genotype_names[snp1_index][0]
            geneindex1_select = genotype_names[snp1_index][1].split(',')
            gene_snp1 = geneindex1_select[0]
            
            job_value_temp["snp1"]=name_snp1   
            job_value_temp["gene1"]=gene_snp1    #drop?????        

            #Check if the snp has been deemed unsuitable (monomorph or not in hwe)
            if len(snps_not_HWE)>0 and name_snp1 == snps_not_HWE[0]:
               snps_not_HWE.pop(0) 
               continue
            for snp2_index in range(snp1_index+1,len(genotype_names)):
#                name_snp2, gene_snp2 = genotype_names[snp2_index]#can also be used
                name_snp2 = genotype_names[snp2_index][0]
                geneindex2_select = genotype_names[snp2_index][1].split(',')
                gene_snp2 = geneindex2_select[0]

                job_value_temp["snp2"]=name_snp2   
                job_value_temp["gene2"]=gene_snp2  #drop?????          

                """Created matrix with selected trait and snps"""
                                
                # construct a 3-column submatrix with a SNP pair and a trait
                data_use_matrix = np.matrix([phenotypes[trait_index], genotypes[snp1_index], genotypes[snp2_index]])
                data_use_matrix = data_use_matrix.flatten().reshape(-1,3, order="F") 
                
                #Removes cases with no trait value or no genotype values for one or both snps
                data_use_matrix = filter_selection(data_use_matrix)

#FILTERING FOR MISSING TRAIT VALUES COULD BE DONE EARLIER, WHICH MAY SAVE SOME COMPUTER TIME. OFTEN THIS MAY BE ONLY A FEW, SO IS IT WORTHWILE, B? HOWEVER IT CAN BE SUBSTANTIAL AND IN RARE INSTANCES NO TRAIT VALUES ARE PRESENT. ALSO, THER MAY BE LESS THAN 10 CASSES LEFT IN WHICH INSTANCE, WHICH IS NOT COMPATIBLE WITH EPISTASIS SIGNIFICANCE STEP! THIS IS ACCOUNTED FOR!

                #In the filtering monomorph can emerge
                #Check if one or both snps have become monomorph

                mono_check = filter_genotypes_select(data_use_matrix,snps_not_HWE,job_value_temp)

                #Return None if not 1 or 2
                if mono_check == 1:
                	print "SNP1 Not OK"#This should only happen for the first monomorph snp
                	break
                if mono_check == 2:
                	print "SNP2 Not OK" #???               	
                	continue

                """Calculate allel frequences and test for HWE"""

                allel_val= get_allelfreqences(data_use_matrix,job_value_temp,cut_off)

#return case_valid, allelfreq, gtall,hsnpA, hsnpB

                """Check for HWE; skip if non-equilibirum"""
                """Select for minimum allele frenquence"""  
                if allel_val[3] == 0:#This is only be relevant for the first snp; if tested before it can be deleted

                   break #OK
                if allel_val[4] == 0 or allel_val[0] == 0:#MAF test
                   snps_not_HWE.append(job_value_temp["snp2"])
                   continue 

                """Calculating means"""
                geno_mean = epistasis_calc_mean(data_use_matrix,allel_val[1])
#return ltraitgt, mtraitgt, vtraitgt, TotalgMeans

                """Calculating significans of epistasis"""
                episign = epivar_sign(geno_mean[0],geno_mean[1],geno_mean[2])
#return Fstat, dfdenom, Fpvalue

                """Check for epistasis significance; skip if non-significant"""
                if episign[2] >= bonf_correct:
                    continue
                
                """Calculating total variance"""
                tot_var = epistasis_calc_totvar(allel_val[1],geno_mean[1],geno_mean[3])
                job_value_temp["PhenVar"] = tot_var[0]
                job_value_temp["TgenVar"] = tot_var[1]
           
#return PhenVar, TGenVar

                """Calculating epistasis"""
                
                epi_calc = epistasis_calc_var(allel_val[1],geno_mean[1],geno_mean[3],episign,tot_var,job_value_temp)
#return gteffects, gtvariance, gtvariance_sumLD,gtvariance_sumNLD

#ENTROPY
#Calculate entropy for each snp using alle_val[2
                entropy_data = entropy_shannon.entropy_shannon_init(allel_val[2])
#return entsumma_all,entsumma_geno,entsumma_geno_exp,data-input as a tuple              
                trait_result["Entropy_data"].append(entropy_data)
#NB this is a list of tuples, whos elements are strings

#Calculate weighted entropy for each snp using alle_val[2] and epi_calc[0]=effect sizes
                entropy_data_weight = entropy_shannon_weight.entropy_shannon_init_weight(allel_val[2],epi_calc[0])
#return entsumma_all as a string      
                trait_result["Entropy_data_weight"].append(entropy_data_weight)
#NB this is a list of strings

                trait_result["snp_interactions"].append((job_value_temp["snp1"],job_value_temp["snp2"], allel_val[2], epi_calc,entropy_data,entropy_data_weight))
                trait_result["AllGenFreq"].append(allel_val[2])
                
#Calculate signifcance
                epi_sum_out.append(epi_calc[2:3][0])  #lostof lists, no header!!LD and NLD

        epistasis_significant = len(epi_sum_out)
        relative_sign = float(epistasis_significant)/float(num_snp_tests)
        trait_result["relative_significant"] = relative_sign 
        trait_result["epistasis_significant"] = epistasis_significant  
        
#ENTROPY SUMMARY        
#Calculate entropy summary values
        entropy_val_sum = entropy_shannon.entropy_shannon_sum(trait_result["Entropy_data"])
        trait_result["Entropy"] = entropy_val_sum

#Calculate weighted entropy summary values
        entropy_val_sum_weight = entropy_shannon_weight.entropy_shannon_sum_w1(trait_result["Entropy_data_weight"])
        trait_result["Entropy_weight"] = entropy_val_sum_weight

        results.append(trait_result)            

#SUMMERING AF RESULTATER FRA GRID. ER VEL ET POST-GRID SERVER PROBLEM?
            
    return results

def entropy_shanonX(allel_data):		
	#Calculate Shanon entropies
	
	#Alleles
	allel_raw = allel_data[2]
	
	# Alleles in først SNP
	allelA1 = float(allel_raw[5])
	allelA2 = 1 - allelA1
	#This should be superflous as monmoprhs has been discarded
	if allelA1 != 0:
		ent_A1 = -allelA1*np.log2(allelA1)
		ent_A2 = -allelA2*np.log2(allelA2)
	else: 
		ent_A1 = 0
		ent_A2 = 0	
	#Sum entropy for snp A alleles
	ent_snpA_all = ent_A1 + ent_A2
	
	# Alleles in second SNP
	allelB1 = float(allel_raw[6])
#	allelB3 = float(allel_raw[7])
	allelB4 = float(allel_raw[8])		

	allelB7 = float(allel_raw[11])	
	allelB8 = float(allel_raw[12])			
	allelB9 = float(allel_raw[13])	
	allelB2 = 1 - allelB1
	#This should be superflous as monmoprhs has been discarded
	if allelB1 != 0:
		ent_B1 = -allelB1*np.log2(allelB1)
		ent_B2 = -allelB2*np.log2(allelB2)
	else: 
		ent_B1 = 0
		ent_B2 = 0
	
	#Sum entropy for snp B alleles
	ent_snpB_all = ent_B1 + ent_B2

	#Epistatic entropy alleles
	ent_epi_allel = ent_snpA_all + ent_snpB_all

	
	#Genotypes counted
	summa_geno = float(allel_raw[2])+float(allel_raw[3])+float(allel_raw[4])
	#Both SNPs must have the same number of genotypes

	# Genotypes in first SNP	
	freq_genoA = [float(allel_raw[2])/summa_geno,float(allel_raw[3])/summa_geno,float(allel_raw[4])/summa_geno]

	if freq_genoA[0] != 0:
		ent_gA11 = -freq_genoA[0]*np.log2(freq_genoA[0])
	else: ent_gA11 = 0
	if freq_genoA[1] != 0:
		ent_gA12 = -freq_genoA[1]*np.log2(freq_genoA[1])
	else: ent_gA12 = 0
	if freq_genoA[2] != 0:
		ent_gA22 = -freq_genoA[2]*np.log2(freq_genoA[2])
	else: ent_gA22 = 0
	
	#Sum entropy for genotypes A
	ent_genoA_all = ent_gA11 + ent_gA12 + ent_gA22
		
	# Genotypes in second SNP	
	freq_genoB = [float(allel_raw[9])/summa_geno,float(allel_raw[10])/summa_geno,float(allel_raw[11])/summa_geno]
	
	if freq_genoB[0] != 0:
		ent_gB11 = -freq_genoB[0]*np.log2(freq_genoB[0])
	else: ent_gB11 = 0
	if freq_genoB[1] != 0:
		ent_gB12 = -freq_genoB[1]*np.log2(freq_genoB[1])
	else: ent_gB12 = 0
	if freq_genoB[2] != 0:
		ent_gB22 = -freq_genoB[2]*np.log2(freq_genoB[2])
	else: ent_gB22 = 0
	
	#Sum entropy for genotypes B
	ent_genoB_all = ent_gB11 + ent_gB12 + ent_gB22

	#Epistatic entropy genotypes counted
	ent_epi_geno = ent_genoA_all + ent_genoB_all

	#Genotypes expected

	exp_genoA = [allelA1*allelA1,2*allelA1*allelA2,allelA2*allelA2]

	#Superflous precaution as allelfrequences cannot be zero
	if exp_genoA[0] != 0:
		eent_gA11 = -exp_genoA[0]*np.log2(exp_genoA[0])
	else: eent_gA11 = 0
	if exp_genoA[1] != 0:
		eent_gA12 = -exp_genoA[1]*np.log2(exp_genoA[1])
	else: eent_gA12 = 0
	if exp_genoA[2] != 0:
		eent_gA22 = -exp_genoA[2]*np.log2(exp_genoA[2])
	else: eent_gA22 = 0
	
	#Sum entropy for genotypes A expected
	eent_genoA_all = eent_gA11 + eent_gA12 + eent_gA22

	exp_genoB = [allelB1*allelB1,2*allelB1*allelB2,allelB2*allelB2]
	#Superflous precaution as allelfrequences cannot be zero
	if exp_genoB[0] != 0:
		eent_gB11 = -exp_genoB[0]*np.log2(exp_genoB[0])
	else: eent_gB11 = 0
	if exp_genoB[1] != 0:
		eent_gB12 = -exp_genoB[1]*np.log2(exp_genoB[1])
	else: eent_gB12 = 0
	if exp_genoB[2] != 0:
		eent_gB22 = -exp_genoB[2]*np.log2(exp_genoB[2])
	else: eent_gB22 = 0
	
	#Sum entropy for genotypes B expected
	eent_genoB_all = eent_gB11 + eent_gB12 + eent_gB22

	#Epistatic entropy genotypes expected
	eent_epi_geno = eent_genoA_all + eent_genoB_all
	
	entsumma_all = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],ent_A1,ent_A2,ent_snpA_all,job_value_temp["snp2"],job_value_temp["gene2"],ent_B1,ent_B2,ent_snpB_all,ent_epi_allel)

	entsumma_geno = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11} {12}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],ent_gA11,ent_gA12,ent_gA22,ent_genoA_all,job_value_temp["snp2"],job_value_temp["gene2"],ent_gB11,ent_gB12,ent_gB22,ent_genoB_all,ent_epi_geno)

	entsumma_geno_exp = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11} {12}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],eent_gA11,eent_gA12,eent_gA22,eent_genoA_all,job_value_temp["snp2"],job_value_temp["gene2"],eent_gB11,eent_gB12,eent_gB22,eent_genoB_all,eent_epi_geno)

	return entsumma_all, entsumma_geno, entsumma_geno_exp

def snp_list_create(data_raw):
	"""Create a list of snps"""
	snp_list_merge = []		
	for epis in range(1,len(data_raw)):
		allel_conv= data_raw[epis].split()
		print type(allel_conv)
		snp_list_merge.append(allel_conv[0])
		snp_list_merge.append(allel_conv[7])
	snp_list = list(set(snp_list_merge))
	snp_list_size = len(snp_list)#to be used in bonferoni
	return snp_list,snp_list_size

#DETTE SCRIPT KOMBINERER VALGTE TRAIT (I RELATION TIL PARTITION) OG DE TO SNPS DVS ER ET GRID-SCRIPT JOB
def filter_selection(selection_matrix): #OK mf

    #Removes cases with no trait value or no genotype values for one or both snps
    flawed_rows = []
    for row in xrange(len(selection_matrix)): #phenotypes can be zero
        if selection_matrix[row,0] < 0 or selection_matrix[row,1] < 1 or selection_matrix[row,2] < 1: 
            flawed_rows.append(row)
    selection_matrix = np.delete(selection_matrix, flawed_rows, 0)
        
    return selection_matrix
    
def get_allelfreqences(data_use_matrix,job_value_temp, cut_off):#OK 

    #sum_cases = data_use_matrix.shape[0]-1
    sum_cases = data_use_matrix.shape[0]# B : column title no longer in first row 

    genotypeA = [0,0,0]
    genotypeB = [0,0,0]
    for gt in xrange(0, len(data_use_matrix)):
        genotypeA[int(data_use_matrix[gt,1])-1] +=1
        genotypeB[int(data_use_matrix[gt,2])-1] +=1

    """Check for monomorph snps. If so, further calculation is skipped"""
  
    
    """Allele frequencies are calculated"""
    allelA = float(2*genotypeA[0] + genotypeA[1])/(sum_cases*2)
    allela = float(2*genotypeA[2] + genotypeA[1])/(sum_cases*2)
    allelB = float(2*genotypeB[0] + genotypeB[1])/(sum_cases*2)
    allelb = float(2*genotypeB[2] + genotypeB[1])/(sum_cases*2)
    allelfreq = list((allelA, allela, allelB, allelb))
    
    """Testing for minimum allel frequence"""
    #First SNP
    minAFA = min(allelA,allela)
    totAFA = allelA+allela
    minA = minAFA/totAFA
    #Second SNP
    minAFB = min(allelB,allelb)
    totAFB = allelB+allelb
    minB = minAFB/totAFB

    case_valid = 1
    if minA<job_value_temp["MAF"] or minB<job_value_temp["MAF"]:
        case_valid = 0    
    """HWE calculations"""
    #Gene A
    obs_sum_geneA = sum(genotypeA)
    exp_geneA = [allelfreq[0]*allelfreq[0]*obs_sum_geneA,
                2*allelfreq[0]*allelfreq[1]*obs_sum_geneA,
                allelfreq[1]*allelfreq[1]*obs_sum_geneA]

    obs_exp_diff_powA = [pow(genotypeA[0]-exp_geneA[0],2),
                         pow(genotypeA[1]-exp_geneA[1],2),
                         pow(genotypeA[2]-exp_geneA[2],2)]
    chis_valA = [obs_exp_diff_powA[0]/exp_geneA[0],
                 obs_exp_diff_powA[1]/exp_geneA[1],
                 obs_exp_diff_powA[2]/exp_geneA[2]]
    chis_val_sumA = sum(chis_valA)
    hwe_geneA = stats.chisqprob(chis_val_sumA,2)

    #Gene B
    obs_sum_geneB = sum(genotypeB)
    exp_geneB = [allelfreq[2]*allelfreq[2]*obs_sum_geneB,
                2*allelfreq[2]*allelfreq[3]*obs_sum_geneB,
                allelfreq[3]*allelfreq[3]*obs_sum_geneB]

    obs_exp_diff_powB = [pow(genotypeB[0]-exp_geneB[0],2),
                         pow(genotypeB[1]-exp_geneB[1],2),
                         pow(genotypeB[2]-exp_geneB[2],2)]
    chis_valB = [obs_exp_diff_powB[0]/exp_geneB[0],
                 obs_exp_diff_powB[1]/exp_geneB[1],
                 obs_exp_diff_powB[2]/exp_geneB[2]]
    chis_val_sumB = sum(chis_valB)
    hwe_geneB = stats.chisqprob(chis_val_sumB,2)
    
    gtall = ""
    hsnpA=1
    hsnpB=1
#HVIS HWE IKKE ER OPFYLDT FOR BEGGE GENER SKIPPES DISSE (EVT MUD UDSKRIFT) OG DER FORTSÆTTES TIL NÆSTE KOMBINATION AF SNPS. GRID-SCRIPT!
#If HWE is not fulfilled..

    if hwe_geneA >= cut_off and hwe_geneB >= cut_off:

        gtall = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11} {12} {13}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],genotypeA[0],
           genotypeA[1],genotypeA[2],
           allelA,allela,job_value_temp["snp2"],job_value_temp["gene2"],genotypeB[0],genotypeB[1],genotypeB[2],
           allelB,allelb)
  #      case_valid = 1
    else:
   #     case_valid = 0
        if hwe_geneA < cut_off:
           hsnpA =0
        else:    
           hsnpA = 1
        if hwe_geneB < cut_off:
           hsnpB =0
        else:    
           hsnpB = 1
     
    return case_valid, allelfreq, gtall, hsnpA, hsnpB


#DETTE ER ET GRID-SCRIPT

#BASALE BEREGNINGER. GRID-SCRIPT

def epistasis_calc_mean(data_use_matrix,allelfreq):

    #Overall mean
    #all_mean = np.average(data_use_matrix[1:,0].astype(float))
    all_mean = np.average(data_use_matrix[:,0].astype(float)) # B : column title no longer i first row
    #To be printed to file containing all the means and variances of 2-snps results

    """Traivalues for every combinations of two snps"""
    traitgt11 = []
    traitgt12 = []
    traitgt13 = []
    traitgt21 = []
    traitgt22 = []
    traitgt23 = []
    traitgt31 = []
    traitgt32 = []
    traitgt33 = []

    """Calculating the trait means for 2-snp interactions"""
    # B : column title no longer in first row
    for gt in range(0, len(data_use_matrix)-1): 
        if data_use_matrix[gt,1].astype(int)==1 and data_use_matrix[gt,2].astype(int)==1:
           traitgt11.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==1 and data_use_matrix[gt,2].astype(int)==2:
           traitgt12.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==1 and data_use_matrix[gt,2].astype(int)==3:
           traitgt13.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==2 and data_use_matrix[gt,2].astype(int)==1:
           traitgt21.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==2 and data_use_matrix[gt,2].astype(int)==2:
           traitgt22.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==2 and data_use_matrix[gt,2].astype(int)==3:
           traitgt23.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==3 and data_use_matrix[gt,2].astype(int)==1:
           traitgt31.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==3 and data_use_matrix[gt,2].astype(int)==2:
           traitgt32.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==3 and data_use_matrix[gt,2].astype(int)==3:
           traitgt33.append(data_use_matrix[gt,0].astype(float))

    """Number of subjects in each group"""
#    global ltraitgt
    ltraitgt = [len(traitgt11),len(traitgt12),len(traitgt13),len(traitgt21),\
                len(traitgt22),len(traitgt23),len(traitgt31),len(traitgt32),\
                len(traitgt33)]
    """Average of traits, trait-menas"""
    mtraitgt11 = np.average(traitgt11)
    mtraitgt12 = np.average(traitgt12)
    mtraitgt13 = np.average(traitgt13)
    mtraitgt21 = np.average(traitgt21)
    mtraitgt22 = np.average(traitgt22)
    mtraitgt23 = np.average(traitgt23)
    mtraitgt31 = np.average(traitgt31)
    mtraitgt32 = np.average(traitgt32)
    mtraitgt33 = np.average(traitgt33)

    mtraitgt = [mtraitgt11,mtraitgt12,mtraitgt13,mtraitgt21,\
                mtraitgt22,mtraitgt23,mtraitgt31,mtraitgt32,mtraitgt33]
    #NB traits without content (0) is nan
    for i in range(len(ltraitgt)):
        if ltraitgt[i] == 0 :mtraitgt[i]=0
#    global vtraitgt
    vtraitgt11 = np.var(traitgt11)
    vtraitgt12 = np.var(traitgt12)
    vtraitgt13 = np.var(traitgt13)
    vtraitgt21 = np.var(traitgt21)
    vtraitgt22 = np.var(traitgt22)
    vtraitgt23 = np.var(traitgt23)
    vtraitgt31 = np.var(traitgt31)
    vtraitgt32 = np.var(traitgt32)
    vtraitgt33 = np.var(traitgt33)
    vtraitgt = [vtraitgt11,vtraitgt12,vtraitgt13,vtraitgt21,\
                vtraitgt22,vtraitgt23,vtraitgt31,vtraitgt32,vtraitgt33]

    """Mean values conditioned on HWE of snps"""

  #Mean values for genotypes calculated as HWE frequences = OK
    #gMean is eq. to gVaacc etc in R-script
    gMean11 = allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[2]*mtraitgt[0]
    gMean12 = 2*allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[3]*mtraitgt[1]
    gMean13 = allelfreq[0]*allelfreq[0]*allelfreq[3]*allelfreq[3]*mtraitgt[2]
    gMean21 = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[3]
    gMean22 = 4*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[4]
    gMean23 = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[5]
    gMean31 = allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[6]
    gMean32 = 2*allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[7]
    gMean33 = allelfreq[1]*allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[8]

    gMean = [gMean11,gMean12,gMean13,gMean21,gMean22,gMean23,gMean31,gMean32,gMean33]
    #TotalgMeans = TotgValues in original R-script
    TotalgMeans = np.sum(gMean)

    return ltraitgt, mtraitgt, vtraitgt, TotalgMeans

def epistasis_calc_totvar(allelfreq,mtraitgt,TotalgMeans):

    #NB change to AABB etc 
    gVar11 = allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[2]*pow(mtraitgt[0],2)
    gVar12 = 2*allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[3]*pow(mtraitgt[1],2)
    gVar13 = allelfreq[0]*allelfreq[0]*allelfreq[3]*allelfreq[3]*pow(mtraitgt[2],2)
    gVar21 = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(mtraitgt[3],2)
    gVar22 = 4*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(mtraitgt[4],2)
    gVar23 = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(mtraitgt[5],2)
    gVar31 = allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(mtraitgt[6],2)
    gVar32 = 2*allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(mtraitgt[7],2)
    gVar33 = allelfreq[1]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(mtraitgt[8],2)

    PhenVar = np.sum([gVar11,gVar12,gVar13,gVar21,gVar22,gVar23,gVar31,gVar32,gVar33])
    #For checksum: TGenVar the same as in R-script
    TGenVar = PhenVar - pow(TotalgMeans,2)

    return PhenVar, TGenVar
    
##NB DETTE ER KERNEN I BEREGNINGERNE OG ALTÅ GRID-SCRIPT


def epistasis_calc_var(allelfreq,mtraitgt,TotalgMeans,episign,tot_var,job_value_temp):

    """Genotypic values"""

#Allele A OK, fix A
    addValAABB = allelfreq[0]*allelfreq[2]*allelfreq[2]*mtraitgt[0]
    addValAABb = 2*allelfreq[0]*allelfreq[2]*allelfreq[3]*mtraitgt[1]
    addValAAbb = allelfreq[0]*allelfreq[3]*allelfreq[3]*mtraitgt[2]
    addValAaBB = allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[3]
    addValAaBb = 2*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[4]
    addValAabb = allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[5]

    addValA = [addValAABB,addValAABb,addValAAbb,addValAaBB,addValAaBb,addValAabb]
    addValATot = np.sum(addValA)
 
#Allel a OK, fix a 
    addValAaBB = allelfreq[0]*allelfreq[2]*allelfreq[2]*mtraitgt[3]
    addValAaBb = 2*allelfreq[0]*allelfreq[2]*allelfreq[3]*mtraitgt[4]
    addValAabb = allelfreq[0]*allelfreq[3]*allelfreq[3]*mtraitgt[5]
    addValaaBB = allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[6]
    addValaaBb = 2*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[7]
    addValaabb = allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[8]

    addVala = [addValAaBB,addValAaBb,addValAabb,addValaaBB,addValaaBb,addValaabb]
    addValaTot = np.sum(addVala)

#Allel B OK, fix B
    addValAABB = allelfreq[0]*allelfreq[0]*allelfreq[2]*mtraitgt[0]
    addValAaBB = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*mtraitgt[3]
    addValaaBB = allelfreq[1]*allelfreq[1]*allelfreq[2]*mtraitgt[6]
    addValAABb = allelfreq[0]*allelfreq[0]*allelfreq[3]*mtraitgt[1]
    addValAaBb = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*mtraitgt[4]
    addValaaBb = allelfreq[1]*allelfreq[1]*allelfreq[3]*mtraitgt[7]

    addValB = [addValAABB,addValAaBB,addValaaBB,addValAABb,addValAaBb,addValaaBb]
    addValBTot = np.sum(addValB)
 
#Allel b OK, fix b
    addValAABb = allelfreq[0]*allelfreq[0]*allelfreq[2]*mtraitgt[1]
    addValAaBb = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*mtraitgt[4]
    addValaaBb = allelfreq[1]*allelfreq[1]*allelfreq[2]*mtraitgt[7]
    addValAAbb = allelfreq[0]*allelfreq[0]*allelfreq[3]*mtraitgt[2]
    addValAabb = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*mtraitgt[5]
    addValaabb = allelfreq[1]*allelfreq[1]*allelfreq[3]*mtraitgt[8]

    addValb = [addValAABb,addValAaBb,addValaaBb,addValAAbb,addValAabb,addValaabb]
    addValbTot = np.sum(addValb)
  
    """Additive effects, zero centered OK"""

    addEffA = addValATot - TotalgMeans
    addEffa = addValaTot - TotalgMeans
    addEffB = addValBTot - TotalgMeans
    addEffb = addValbTot - TotalgMeans
    #Could this be used to say something about increasing/decreasing the trait?
    addCheckSNP1 = [addEffA,addEffa]
    addCheckSNP2 = [addEffB,addEffb]

    """Additive variances"""
    AddVarSNP1 = 2*(allelfreq[0]*pow(addEffA,2) + allelfreq[1]*pow(addEffa,2))
    AddVarSNP2 = 2*(allelfreq[2]*pow(addEffB,2) + allelfreq[3]*pow(addEffb,2))

    """Total additive variance OK"""
    AddVarTot = AddVarSNP1 + AddVarSNP2

    """Dominant values, effects and variance"""

#Conditonal genotypic values SNP1 OK
    domValAA = pow(allelfreq[2],2)*mtraitgt[0] + 2*allelfreq[2]*allelfreq[3]*mtraitgt[1] + pow(allelfreq[3],2)*mtraitgt[2]
    domValAa = pow(allelfreq[2],2)*mtraitgt[3] + 2*allelfreq[2]*allelfreq[3]*mtraitgt[4] + pow(allelfreq[3],2)*mtraitgt[5]
    domValaa = pow(allelfreq[2],2)*mtraitgt[6] + 2*allelfreq[2]*allelfreq[3]*mtraitgt[7] + pow(allelfreq[3],2)*mtraitgt[8]
#Conditonal genotypic values SNP2 OK
    domValBB = pow(allelfreq[0],2)*mtraitgt[0] + 2*allelfreq[0]*allelfreq[1]*mtraitgt[3] + pow(allelfreq[1],2)*mtraitgt[6]
    domValBb = pow(allelfreq[0],2)*mtraitgt[1] + 2*allelfreq[0]*allelfreq[1]*mtraitgt[4] + pow(allelfreq[1],2)*mtraitgt[7]
    domValbb = pow(allelfreq[0],2)*mtraitgt[2] + 2*allelfreq[0]*allelfreq[1]*mtraitgt[5] + pow(allelfreq[1],2)*mtraitgt[8]

#Dominant effects OK
    domEffAA = domValAA - TotalgMeans - addEffA - addEffA
    domEffAa = domValAa - TotalgMeans - addEffA - addEffa
    domEffaa = domValaa - TotalgMeans - addEffa - addEffa
    domEffBB = domValBB - TotalgMeans - addEffB - addEffB
    domEffBb = domValBb - TotalgMeans - addEffB - addEffb
    domEffbb = domValbb - TotalgMeans - addEffb - addEffb
    domCheckSNP1 = [domEffAA,domEffAa,domEffaa]
    domCheckSNP2 = [domEffBB,domEffBb,domEffbb]
    
#Dominant variance weigthed OK
    domVarAA = allelfreq[0]*allelfreq[0]*pow(domEffAA,2)
    domVarAa = 2*allelfreq[0]*allelfreq[1]*pow(domEffAa,2)
    domVaraa = allelfreq[1]*allelfreq[1]*pow(domEffaa,2)
    domVarBB = allelfreq[2]*allelfreq[2]*pow(domEffBB,2)
    domVarBb = 2*allelfreq[2]*allelfreq[3]*pow(domEffBb,2)
    domVarbb = allelfreq[3]*allelfreq[3]*pow(domEffbb,2)
    vardomCheckSNP1 = [domVarAA,domVarAa,domVaraa]
    vardomCheckSNP2 = [domVarBB,domVarBb,domVarbb]

    domVarTot = np.sum(vardomCheckSNP1 + vardomCheckSNP2)

    """Additive*additive values, effects and variance"""

#Additive*Additive values OK

    addaddValAB = allelfreq[0]*allelfreq[2]*mtraitgt[0]+allelfreq[0]*allelfreq[3]*mtraitgt[1]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[3]+allelfreq[1]*allelfreq[3]*mtraitgt[4]
    addaddValAb = allelfreq[0]*allelfreq[2]*mtraitgt[1]+allelfreq[0]*allelfreq[3]*mtraitgt[2]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[4]+allelfreq[1]*allelfreq[3]*mtraitgt[5]
    addaddValaB = allelfreq[0]*allelfreq[2]*mtraitgt[3]+allelfreq[0]*allelfreq[3]*mtraitgt[4]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[6]+allelfreq[1]*allelfreq[3]*mtraitgt[7]
    addaddValab = allelfreq[0]*allelfreq[2]*mtraitgt[4]+allelfreq[0]*allelfreq[3]*mtraitgt[5]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[7]+allelfreq[1]*allelfreq[3]*mtraitgt[8]

    addaddCheckGene = [addaddValAB,addaddValAb,addaddValaB,addaddValab]

#Additive*Additive effects OK
    addaddEffAB = addaddValAB - TotalgMeans - addEffA - addEffB
    addaddEffAb = addaddValAb - TotalgMeans - addEffA - addEffb
    addaddEffaB = addaddValaB - TotalgMeans - addEffa - addEffB
    addaddEffab = addaddValab - TotalgMeans - addEffa - addEffb
    
    addaddEff = [addaddEffAB,addaddEffAb,addaddEffaB,addaddEffab]
 
#Additive*Additive variance OK

    addaddVarAB = allelfreq[0]*allelfreq[2]*pow(addaddEffAB,2)
    addaddVarAb = allelfreq[0]*allelfreq[3]*pow(addaddEffAb,2)
    addaddVaraB = allelfreq[1]*allelfreq[2]*pow(addaddEffaB,2)
    addaddVarab = allelfreq[1]*allelfreq[3]*pow(addaddEffab,2)

    addaddVar = [addaddVarAB,addaddVarAb,addaddVaraB,addaddVarab]
    addaddVarTot = 4*np.sum(addaddVar)

    """Additive*dominance values, effects and variances"""

#Additive*dominance values OK
    #SNP1
    adddomValAxBB = allelfreq[0]*mtraitgt[0] + allelfreq[1]*mtraitgt[3]
    adddomValAxBb = allelfreq[0]*mtraitgt[1] + allelfreq[1]*mtraitgt[4]
    adddomValAxbb = allelfreq[0]*mtraitgt[2] + allelfreq[1]*mtraitgt[5]
    adddomValaxBB = allelfreq[0]*mtraitgt[3] + allelfreq[1]*mtraitgt[6]
    adddomValaxBb = allelfreq[0]*mtraitgt[4] + allelfreq[1]*mtraitgt[7]
    adddomValaxbb = allelfreq[0]*mtraitgt[5] + allelfreq[1]*mtraitgt[8]

    adddomValA = [adddomValAxBB,adddomValAxBb,adddomValAxbb,\
                  adddomValaxBB,adddomValaxBb,adddomValaxbb]
 
    #SNP2
    adddomValAABx = allelfreq[2]*mtraitgt[0] + allelfreq[3]*mtraitgt[1]
    adddomValAaBx = allelfreq[2]*mtraitgt[3] + allelfreq[3]*mtraitgt[4]    
    adddomValaaBx = allelfreq[2]*mtraitgt[6] + allelfreq[3]*mtraitgt[7]
    adddomValAAbx = allelfreq[2]*mtraitgt[1] + allelfreq[3]*mtraitgt[2]
    adddomValAabx = allelfreq[2]*mtraitgt[4] + allelfreq[3]*mtraitgt[5]
    adddomValaabx = allelfreq[2]*mtraitgt[7] + allelfreq[3]*mtraitgt[8]

    adddomValB = [adddomValAABx,adddomValAaBx,adddomValaaBx,\
                  adddomValAAbx,adddomValAabx,adddomValaabx]
    
#Additive*dominance effects OK
    #SNP1
    adddomEffAxBB = adddomValAxBB - TotalgMeans - addEffA - addEffB - addEffB\
                    -domEffBB - addaddEffAB - addaddEffAB
    adddomEffAxBb = adddomValAxBb - TotalgMeans - addEffA - addEffB - addEffb\
                    -domEffBb - addaddEffAB - addaddEffAb
    adddomEffAxbb = adddomValAxbb - TotalgMeans - addEffA - addEffb - addEffb\
                    -domEffbb - addaddEffAb - addaddEffAb
    adddomEffaxBB = adddomValaxBB - TotalgMeans - addEffa - addEffB - addEffB\
                    -domEffBB - addaddEffaB - addaddEffaB
    adddomEffaxBb = adddomValaxBb - TotalgMeans - addEffa - addEffB - addEffb\
                    -domEffBb - addaddEffaB - addaddEffab
    adddomEffaxbb = adddomValaxbb - TotalgMeans - addEffa - addEffb - addEffb\
                    -domEffbb - addaddEffab - addaddEffab

    adddomEffA = [adddomEffAxBB,adddomEffAxBb,adddomEffAxbb,adddomEffaxBB,adddomEffaxBb,adddomEffaxbb]
 
    #SNP2
    adddomEffAABx = adddomValAABx - TotalgMeans - addEffA - addEffA - addEffB\
                    -domEffAA - addaddEffAB - addaddEffAB
    adddomEffAaBx = adddomValAaBx - TotalgMeans - addEffA - addEffa - addEffB\
                    -domEffAa - addaddEffAB - addaddEffaB
    adddomEffaaBx = adddomValaaBx - TotalgMeans - addEffa - addEffa - addEffB\
                    -domEffaa - addaddEffaB - addaddEffaB
    adddomEffAAbx = adddomValAAbx - TotalgMeans - addEffA - addEffA - addEffb\
                    -domEffAA - addaddEffAb - addaddEffAb
    adddomEffAabx = adddomValAabx - TotalgMeans - addEffA - addEffa - addEffb\
                    -domEffAa - addaddEffAb - addaddEffab
    adddomEffaabx = adddomValaabx - TotalgMeans - addEffa - addEffa - addEffb\
                    -domEffaa - addaddEffab - addaddEffab

    adddomEffB = [adddomEffAABx,adddomEffAaBx,adddomEffaaBx,adddomEffAAbx,adddomEffAabx,adddomEffaabx]

#Additive*dominance variance OK
    #SNP1
    addomVarAxBB = allelfreq[0]*allelfreq[2]*allelfreq[2]*pow(adddomEffAxBB,2)
    addomVarAxBb = 2*allelfreq[0]*allelfreq[2]*allelfreq[3]*pow(adddomEffAxBb,2)
    addomVarAxbb = allelfreq[0]*allelfreq[3]*allelfreq[3]*pow(adddomEffAxbb,2)
    addomVaraxBB = allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(adddomEffaxBB,2)
    addomVaraxBb = 2*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(adddomEffaxBb,2)
    addomVaraxbb = allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(adddomEffaxbb,2)

    adddomVarA = [addomVarAxBB,addomVarAxBb,addomVarAxbb,addomVaraxBB,addomVaraxBb,addomVaraxbb]
 
    #SNP2
    addomVarAABx = allelfreq[0]*allelfreq[0]*allelfreq[2]*pow(adddomEffAABx,2)
    addomVarAaBx = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*pow(adddomEffAaBx,2)
    addomVaraaBx = allelfreq[1]*allelfreq[1]*allelfreq[2]*pow(adddomEffaaBx,2)
    addomVarAAbx = allelfreq[0]*allelfreq[0]*allelfreq[3]*pow(adddomEffAAbx,2)
    addomVarAabx = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*pow(adddomEffAabx,2)
    addomVaraabx = allelfreq[1]*allelfreq[1]*allelfreq[3]*pow(adddomEffaabx,2)
        
    adddomVarB = [addomVarAABx,addomVarAaBx,addomVaraaBx,addomVarAAbx,addomVarAabx,addomVaraabx]

    adddomVarTotA = 2*np.sum(adddomVarA)
    adddomVarTotB = 2*np.sum(adddomVarB)
    adddomVarTot = adddomVarTotA+adddomVarTotB 
    
##Dominant*dominance effects and variances (no intermediary values)

#Dominance*dominance effects OK

    domdomEffAABB = mtraitgt[0]-TotalgMeans-addEffA-addEffA-addEffB-addEffB\
            -domEffAA-domEffBB-addaddEffAB-addaddEffAB-addaddEffAB-addaddEffAB\
            -adddomEffAxBB-adddomEffAxBB-adddomEffAABx-adddomEffAABx
    domdomEffAABb = mtraitgt[1]-TotalgMeans-addEffA-addEffA-addEffB-addEffb\
            -domEffAA-domEffBb-addaddEffAB-addaddEffAB-addaddEffAb-addaddEffAb\
            -adddomEffAxBb-adddomEffAxBb-adddomEffAABx-adddomEffAAbx
    domdomEffAAbb = mtraitgt[2]-TotalgMeans-addEffA-addEffA-addEffb-addEffb\
            -domEffAA-domEffbb-addaddEffAb-addaddEffAb-addaddEffAb-addaddEffAb\
            -adddomEffAxbb-adddomEffAxbb-adddomEffAAbx-adddomEffAAbx
    domdomEffAaBB = mtraitgt[3]-TotalgMeans-addEffA-addEffa-addEffB-addEffB\
            -domEffAa-domEffBB-addaddEffAB-addaddEffAB-addaddEffaB-addaddEffaB\
            -adddomEffAxBB-adddomEffaxBB-adddomEffAaBx-adddomEffAaBx
    domdomEffAaBb = mtraitgt[4]-TotalgMeans-addEffA-addEffa-addEffB-addEffb\
            -domEffAa-domEffBb-addaddEffAB-addaddEffAb-addaddEffaB-addaddEffab\
            -adddomEffAxBb-adddomEffaxBb-adddomEffAaBx-adddomEffAabx
    domdomEffAabb = mtraitgt[5]-TotalgMeans-addEffA-addEffa-addEffb-addEffb\
            -domEffAa-domEffbb-addaddEffAb-addaddEffAb-addaddEffab-addaddEffab\
            -adddomEffAxbb-adddomEffaxbb-adddomEffAabx-adddomEffAabx
    domdomEffaaBB = mtraitgt[6]-TotalgMeans-addEffa-addEffa-addEffB-addEffB\
            -domEffaa-domEffBB-addaddEffaB-addaddEffaB-addaddEffaB-addaddEffaB\
            -adddomEffaxBB-adddomEffaxBB-adddomEffaaBx-adddomEffaaBx
    domdomEffaaBb = mtraitgt[7]-TotalgMeans-addEffa-addEffa-addEffB-addEffb\
            -domEffaa-domEffBb-addaddEffaB-addaddEffaB-addaddEffab-addaddEffab\
            -adddomEffaxBb-adddomEffaxBb-adddomEffaaBx-adddomEffaabx
    domdomEffaabb = mtraitgt[8]-TotalgMeans-addEffa-addEffa-addEffb-addEffb\
            -domEffaa-domEffbb-addaddEffab-addaddEffab-addaddEffab-addaddEffab\
            -adddomEffaxbb-adddomEffaxbb-adddomEffaabx-adddomEffaabx


    domdomEff = [domdomEffAABB,domdomEffAABb,domdomEffAAbb,domdomEffAaBB,\
                 domdomEffAaBb,domdomEffAabb,domdomEffaaBB,domdomEffaaBb,domdomEffaabb]
    
#Dominance*dominance variance OK
    domdomVarAABB = allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[2]*pow(domdomEffAABB,2)
    domdomVarAABb = 2*allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[3]*pow(domdomEffAABb,2)
    domdomVarAAbb = allelfreq[0]*allelfreq[0]*allelfreq[3]*allelfreq[3]*pow(domdomEffAAbb,2)
    domdomVarAaBB = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(domdomEffAaBB,2)
    domdomVarAaBb = 4*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(domdomEffAaBb,2)
    domdomVarAabb = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(domdomEffAabb,2)
    domdomVaraaBB = allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(domdomEffaaBB,2)
    domdomVaraaBb = 2*allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(domdomEffaaBb,2)
    domdomVaraabb = allelfreq[1]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(domdomEffaabb,2)

    
    domdomVar = [domdomVarAABB,domdomVarAABb,domdomVarAAbb,domdomVarAaBB,\
                 domdomVarAaBb,domdomVarAabb,domdomVaraaBB,domdomVaraaBb,domdomVaraabb]
    domdomVarTot = np.sum(domdomVar)

    gteffects = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11}  \
{12} {13} {14} {15} {16} {17} {18} {19} {20} {21} {22} {23} {24} {25} \
{26} {27} {28} {29} {30} {31} {32} {33} {34} {35} {36} {37} {38}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],addEffA,
           addEffa,domEffAA,domEffAa,domEffaa,job_value_temp["snp2"],job_value_temp["gene2"],addEffB,addEffb,
           domEffBB,domEffBb,domEffbb,addaddEffAB,addaddEffAb,addaddEffaB,addaddEffab,
           adddomEffA[0],adddomEffA[1],adddomEffA[2],adddomEffA[3],adddomEffA[4],adddomEffA[5],
           adddomEffB[0],adddomEffB[1],adddomEffB[2],adddomEffB[3],adddomEffB[4],adddomEffB[5],
           domdomEff[0],domdomEff[1],domdomEff[2],domdomEff[3],domdomEff[4],
           domdomEff[5],domdomEff[6],domdomEff[7],domdomEff[8])
    
    TotGenVar = np.sum(AddVarTot+domVarTot+addaddVarTot+adddomVarTot+domdomVarTot)

    CheckGenVar = tot_var[1] - TotGenVar#kan droppe når konfidence

    gtvariance = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9}  {10} {11} {12} {13} \
{14} {15} {16} {17} {18} {19} {20} {21} {22} {23} {24} {25} \
{26} {27} {28} {29} {30} {31} {32} {33} {34} {35} {36}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],
           AddVarSNP1,domVarAA,domVarAa,domVaraa,job_value_temp["snp2"],job_value_temp["gene2"],AddVarSNP2,
           domVarBB,domVarBb,domVarbb,addaddVarAB,addaddVarAb,addaddVaraB,
           addaddVarab,addomVarAxBB,addomVarAxBb,addomVarAxbb,addomVaraxBB,
           addomVaraxBb,addomVaraxbb,addomVarAABx,addomVarAaBx,addomVaraaBx,
           addomVarAAbx,addomVarAabx,addomVaraabx,domdomVarAABB,domdomVarAABb,
           domdomVarAAbb,domdomVarAaBB,domdomVarAaBb,domdomVarAabb,domdomVaraaBB,
           domdomVaraaBb,domdomVaraabb)
     

 #Check can be deleted when we are confident!!   

    """Summary of variances"""
    epivar = np.sum(addaddVarTot+adddomVarTotA+adddomVarTotB+domdomVarTot)
    relGP = job_value_temp["TgenVar"]/job_value_temp["PhenVar"]
    relEP = epivar/job_value_temp["PhenVar"]
    relEG = epivar/job_value_temp["TgenVar"]


    gtvariance_sumLD = ""
    gtvariance_sumNLD = ""    

    #Only SNPs in different genes.

 #   if job_value_temp["gene1"] != job_value_temp["gene2"] or job_value_temp["gene1"]== "EXOGENIC" or job_value_temp["gene2"]== "EXOGENIC" or job_value_temp["gene1"]== "NA" or job_value_temp["gene2"]== "NA" :

    if job_value_temp["gene1"] != job_value_temp["gene2"]:

        gtvariance_sumNLD = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10}  {11} {12} {13} \
{14} {15} {16} {17} {18} {19}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],job_value_temp["snp2"],job_value_temp["gene2"],job_value_temp["PhenVar"],job_value_temp["TgenVar"],       AddVarTot,domVarTot,addaddVarTot,adddomVarTotA,adddomVarTotB,                                   
           domdomVarTot,epivar,relGP,relEP,relEG,CheckGenVar,episign[0],
           episign[1],episign[2])

    #Only SNPs in the same gene (LD)
    elif job_value_temp["gene1"] == job_value_temp["gene2"]:
    
        gtvariance_sumLD = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10}  {11} {12} {13} \
{14} {15} {16} {17} {18} {19}\n".format(job_value_temp["snp1"],job_value_temp["gene1"],job_value_temp["snp2"],job_value_temp["gene2"],job_value_temp["PhenVar"],job_value_temp["TgenVar"],       AddVarTot,domVarTot,addaddVarTot,adddomVarTotA,adddomVarTotB,                                   
           domdomVarTot,epivar,relGP,relEP,relEG,CheckGenVar,episign[0],
           episign[1],episign[2])

    return gteffects, gtvariance, gtvariance_sumLD,gtvariance_sumNLD


#BASALT GRID-SCRIPT

def epivar_sign(ltraitgt,mtraitgt,vtraitgt):
#significance tests
    #mtraitgt mean of groups NBNB 0-værdier indregnes!! dvs average falder!
    """Average of simple genotypic means"""
    gmaverage = np.average(mtraitgt)
#    print gmaverage
    """Calculation of non-epistatic values"""

    gAA = np.average([mtraitgt[0],mtraitgt[1],mtraitgt[2]])
    gAa = np.average([mtraitgt[3],mtraitgt[4],mtraitgt[5]])
    gaa = np.average([mtraitgt[6],mtraitgt[7],mtraitgt[8]])
    gBB = np.average([mtraitgt[0],mtraitgt[3],mtraitgt[6]])
    gBb = np.average([mtraitgt[1],mtraitgt[4],mtraitgt[7]])
    gbb = np.average([mtraitgt[2],mtraitgt[5],mtraitgt[8]])

    """Non-epistatic genotypic values """

    neAABB = gAA + gBB -  gmaverage
    neAABb = gAA + gBb -  gmaverage
    neAAbb = gAA + gbb -  gmaverage
    neAaBB = gAa + gBB -  gmaverage
    neAaBb = gAa + gBb -  gmaverage
    neAabb = gAa + gbb -  gmaverage
    neaaBB = gaa + gBB -  gmaverage
    neaaBb = gaa + gBb -  gmaverage
    neaabb = gaa + gbb -  gmaverage
    nevect = [neAABB,neAABb,neAAbb,neAaBB,neAaBb,neAabb,neaaBB,neaaBb,neaabb]

    """Calculating significance"""

    evect = []
    varmean =[]
    for i in range(0,8):
        eve = ltraitgt[i]*pow(mtraitgt[i]- nevect[i],2)
        evect.append(eve)
        varm = ltraitgt[i]*vtraitgt[i]
        varmean.append(varm)

    EMS = np.sum(evect)/4
    RMS = np.sum(varmean)/sum(ltraitgt)
    dfdenom = sum(ltraitgt)-9
    if dfdenom <= 0:
       Fstat = "Less than 10 cases"
       dfdenom = 0
       Fpvalue = 0
       return Fstat, dfdenom, Fpvalue
    if RMS == 0:#only if var is zero in all groups, which is a rare, but a possible event
       Fstat = "F-statistics cannot be calculated"
       dfdenom = 0
       Fpvalue = 0
       return Fstat, dfdenom, Fpvalue

    Fstat = EMS/RMS
    """F-stat significance = stats.fprob(4,dfdenom,Fstat)
       From scipy.special.fdtrc is (df nomo,df denom, F-value)"""
    Fpvalue = stats.fprob(4,dfdenom,Fstat)
    return Fstat, dfdenom, Fpvalue

   
#GARBAGE 
#def get_epi_datarawc():
 #   epirawfile = "data/M101 comp  140210W.csv"
  #  xxx=csv.reader(open(epirawfile,"rb"))
   # print(epirawfile)
    #print(xxx)
  #  type(xxx)
   # csvfile = open(epirawfile, "rb").read()

    #dialect = csv.Sniffer().sniff(epirawfile)
  #  print(csvfile)
   # print("number of items 1")
    
